

* Valundersökningen (Swedish Election Study)
/*
List of matched variables
gen decrease_public_sector = VU10_V571			outcome1
gen decrease_defense_expenditure = VU10_V572	outcome2
gen lower_taxes = VU10_V574						outcome4
gen sell_public_enterprise = VU10_V577			outcome6
gen decrease_inequality = VU10_V578				outcome7
gen increase_private_healthcare = VU10_V579		outcome8
gen increase_support_countryside = VU10_V580	outcome13
gen abolish_nuclear = VU10_V581					outcome19
gen leave_eu = VU10_V835						outcome30
gen instate_euro = VU10_V836					outcome31
gen join_nato = VU10_V837						outcome32
gen free_trade = VU10_V838						outcome33
gen decrease_aid = VU10_V839					outcome25
gen less_refugees = VU10_V840					outcome26
gen support_immigrant_culture = VU10_V841		outcome27
gen more_skilled_immigration = VU10_V842		outcome23
gen no_afghan_mission = VU10_V843				outcome34 (support war on terror)
gen six_hour_workday = VU10_V865				outcome14
gen ban_pornography = VU10_V866					outcome15
gen harder_punishment = VU10_V868				outcome17
gen lower_maxtaxa = VU10_V874					outcome10 (lower maxtaxa)
*/



clear all

capture program drop store_est
program define store_est
	local indepvar = "`1'"
	matrix mat_`indepvar'_n[`2',`3']=e(N)
	matrix mat_`indepvar'_beta[`2',`3']=_b[temp_x]
	matrix mat_`indepvar'_se[`2',`3']=_se[temp_x]
	matrix mat_`indepvar'_var[`2',`3']=_se[temp_x]^2
	matrix mat_`indepvar'_t[`2',`3']=_b[temp_x]/_se[temp_x]
end

* SE of effect estimate average
capture program drop calc_SE
program define calc_SE

	local indepvar = "`1'"
	local model = `2'

	* Total sum of variances
	egen sum_var_`indepvar'_`model'=sum(mat_`indepvar'_var`model')
	
	* Second term
	local tempsum = 0
	local actualn = 0
	forvalues i=1(1)`3' {
		if mat_`indepvar'_se[`i',`model']!=. {
			local actualn = `actualn'+1
		}
		forvalues j=1(1)`3' {
			if `i'!=`j' & mat_`indepvar'_se[`i',`model']!=. & mat_`indepvar'_se[`j',`model']!=. {
				local tempsum = `tempsum' + `4'corr_mat[`i',`j'] * mat_`indepvar'_se[`i',`model'] * ///
					mat_`indepvar'_se[`j',`model'] * `4'flip_`indepvar'[`i']  * `4'flip_`indepvar'[`j']
			}
		}
	}
	
	* Total standard error
	gen se_`indepvar'_`model' = sqrt(`tempsum' + sum_var_`indepvar'_`model'[1])*(1/`actualn') in 1

end

cd "C:\Userdata\Shared\Dofiles\PrelDoFiles\PrelRafael\Quantifying bias\replication_code"
use "C:\Userdata\Shared\Dofiles\PrelDoFiles\PrelRafael\Quantifying bias\ext_data\SND0876_VU2010.DTA" 

* Gen dependents (total: 21 included)
gen decrease_public_sector = VU10_V571
gen decrease_defense_expenditure = VU10_V572
gen lower_taxes = VU10_V574
gen sell_public_enterprise = VU10_V577
gen decrease_inequality = VU10_V578
gen increase_private_healthcare = VU10_V579
gen increase_support_countryside = VU10_V580
gen abolish_nuclear = VU10_V581
gen leave_eu = VU10_V835
gen instate_euro = VU10_V836
gen join_nato = VU10_V837
gen free_trade = VU10_V838
gen decrease_aid = VU10_V839
gen less_refugees = VU10_V840
gen support_immigrant_culture = VU10_V841
gen more_skilled_immigration = VU10_V842
gen no_afghan_mission = VU10_V843
gen six_hour_workday = VU10_V865
gen ban_pornography = VU10_V866
gen harder_punishment = VU10_V868
gen lower_maxtaxa = VU10_V874

local j=1
foreach i of varlist decrease_public_sector - lower_maxtaxa {
	gen outcome`j' = `i'
	label variable outcome`j' "`i'"
	local j=`j'+1
	drop `i'
}


* Gen independents
gen Sun2000niva = VU10_V7046
gen educationYears=.
replace educationYears=7 if Sun2000niva==1
replace educationYears=9 if Sun2000niva==2
replace educationYears=12 if Sun2000niva==3
replace educationYears=13 if Sun2000niva==4
replace educationYears=15 if Sun2000niva==5
replace educationYears=19 if Sun2000niva==6
gen bin_edu=.
replace bin_edu=0 if educationYears<=12
replace bin_edu=1 if educationYears>12&educationYears!=.
gen workIncome = VU10_V7045

* Gen controls
gen sex = VU10_S1
gen age_cat = VU10_V7032
gen municipality = VU10_S71
gen rural_city = VU10_V7034
gen occupation = VU10_V7039
gen sector = VU10_V7044

* Inclusion
gen include_indicator=.
replace include_indicator=1 if sex!=. & age_cat!=. & municipality!=. & occupation !=.

* Cleanup
keep outcome1 - include_indicator


* Create necessary matrices for storage
foreach indepvar in educationYears bin_edu workIncome {
	matrix mat_`indepvar'_beta=J(22,2,.)
	matrix mat_`indepvar'_n=J(22,2,.)
	matrix mat_`indepvar'_se=J(22,2,.)
	matrix mat_`indepvar'_var=J(22,2,.)
	matrix mat_`indepvar'_t=J(22,2,.)
}

* Correlation matrix
pwcorr outcome1-outcome21 if include_indicator==1
matrix corr_mat=r(C)


* Run analyses
local i=1
foreach depvar of varlist outcome1-outcome21 {
	
	foreach indepvar in educationYears bin_edu workIncome {
	
		di `i' ", `depvar' `indepvar'"
	
		qui:su `depvar' if `indepvar'!=. & include_indicator!=.
		qui:gen temp_y = (`depvar'-r(mean))/r(sd)
		qui:su `indepvar' if `depvar'!=. & include_indicator!=.
		qui:gen temp_x = (`indepvar'-r(mean))/r(sd)
	
		* Empty model
		qui: reghdfe temp_y temp_x c.sex##i.age_cat if include_indicator==1, noabsorb
		store_est `indepvar' `i' 1
		
		* Naive model
		if "`indepvar'" == "bin_edu"|"`indepvar'" == "educationYears" {
			qui: reghdfe temp_y temp_x c.sex##i.age_cat i.occupation workIncome if include_indicator==1, absorb(municipality) cluster(municipality)
			store_est `indepvar' `i' 2
		}
		else {
			qui: reghdfe temp_y temp_x c.sex##i.age_cat i.occupation educationYears if include_indicator==1, absorb(municipality) cluster(municipality)
			store_est `indepvar' `i' 2
		}
		
		drop temp_x temp_y
		
	}
	local i=`i'+1
}

drop _all
foreach indepvar in educationYears bin_edu workIncome {

	* Turn matrices into variables
	svmat mat_`indepvar'_n
	svmat mat_`indepvar'_beta
	svmat mat_`indepvar'_t
	svmat mat_`indepvar'_se
	svmat mat_`indepvar'_var
	
	* Flip vector
	gen flip_`indepvar' = 1
	replace flip_`indepvar' = -1 if mat_`indepvar'_beta1<0
	
	* Flipped betas
	gen `indepvar'_1 = mat_`indepvar'_beta1 * flip_`indepvar'
	gen `indepvar'_2 = mat_`indepvar'_beta2 * flip_`indepvar'
	
	* Means
	egen m_`indepvar'_1 = mean(`indepvar'_1)
	egen m_`indepvar'_2 = mean(`indepvar'_2)
	gen `indepvar'_m =.
	replace `indepvar'_m = m_`indepvar'_1[1] in 1
	replace `indepvar'_m = m_`indepvar'_2[1] in 2
	drop m_*
	
	* Standard errors
	calc_SE `indepvar' 1 21
	calc_SE `indepvar' 2 21
	gen `indepvar'_se =.
	replace `indepvar'_se = se_`indepvar'_1[1] in 1
	replace `indepvar'_se = se_`indepvar'_2[1] in 2

	gen `indepvar'_Nmin = .
	gen `indepvar'_Nmax = .
	gen `indepvar'_delta = .
	gen `indepvar'_tdelta = .
	forvalues i = 1(1)2 {
		egen Nmin_`indepvar'_`i' = min(mat_`indepvar'_n`i')
		egen Nmax_`indepvar'_`i' = max(mat_`indepvar'_n`i')
		replace `indepvar'_Nmin = Nmin_`indepvar'_`i'[1] in `i'
		replace `indepvar'_Nmax = Nmax_`indepvar'_`i'[1] in `i'
		drop Nmin_`indepvar'_`i' Nmax_`indepvar'_`i'
	}
	replace `indepvar'_delta = `indepvar'_m[2] - `indepvar'_m[1] in 2
	replace `indepvar'_tdelta = `indepvar'_delta[2]/sqrt(`indepvar'_se[1]^2 + `indepvar'_se[2]^2) in 2
	
}

keep *_m *_se *_Nmin *Nmax *_delta *_tdelta
keep if [_n]<3
rename * VU_*
save VU_est.dta, replace


